NovaSAR-1 CEOS-ARD RTC Gamma0 data No description has been provided for this image

Index¶

  • Overview
  • Setup (imports, defaults, dask, odc)
  • Example query
  • Product definition
  • Quality layer
  • Create and apply a good quality pixel mask
  • Plot and browse the data

Overview¶

This notebook demonstrates how to load and use NovaSAR-1 Radiometric Terrain Corrected (RTC) Gamma0 data from the CSIRO NovaSAR Facility.

These CEOS-ARD NovaSAR gamma-0 backscatter data are processed from NovaSAR-1 scenes using the GAMMA software.

Using NovaSAR-1 backscatter data¶

An excellent introduction and overview to using SAR data is provided in the CEOS Laymans SAR Interpretation Guide. This guide has also been converted to a set of set of Jupyter notebooks that you can download from https://github.com/AMA-Labs/cal-notebooks/tree/main/examples/SAR.

Synthetic Aperture Radar operates in the microwave range of the electromagnetic spectrum as an active pulse sent by the satellite and scattered by features on the Earth's surface. The return signal from the surface is measured at the satellite in terms of the signal intensity, phase and polarisation compared to the signal that was sent.

The SAR instrument on the NovaSAR-1 satellite operate in the S-band at approximately 9.4 cm wavelength. This means that it can "see" objects of about this size and larger, and smaller objects are relatively transparent. This makes Sentinel-1 more sensitive to vegetation structural types, sparse and low biomass vegetation, and surface water (smooth and wind affected). Compared to Sentinel-1 (C-band, 5.6 cm), NovaSAR S-band is less sensitive to tree canopies and has better ground penetration.

The SAR signal responds to the orientation and scattering from surface features of comparable size or larger than the wavelength.

  • A bright backscatter value typically means the surface was orientated perpendicular to the signal incidence angle and most of the signal was reflected back to the satellite (direct backscatter)
  • A dark backscatter value means most of the signal was reflected away from the satellite (forward scattering) and typically responds to a smooth surface (relative to the wavelength) such as calm water or bare soil
  • Rough surfaces (relative to the wavelength) result in diffuse scattering where some of the signal is returned to the satellite.
  • Complex surfaces may result in volume scattering (scattering within a tree canopy) or double-bounce scattering (perpendicular objects such as buildings and structures)
  • The relative backscatter values of co-polarisation (HH, VV) and cross-polarisation (HV) measurements can provide information on the scattering characteristics of the surface features.

Using NovaSAR-1 backscatter data requires interpretation of the data for different surface features, including as these features change spatially or in time. It may also be necessary to carefully consider the incidence angle of the SAR signal relative to the surface features using the incidence_angle band or the satellite direction metadata (descending = north to south; ascending = south to north).

Units and conversions¶

The novasar_ceosard_* data are given in Amplitude units. Amplitude can be converted to decibel (dB) or intensity, and vice-versa, with the following equations. Practical Xarray examples are given below.

Amplitude to/from dB:

       dB = 20 * log10(intensity) + K
intensity = 10^((dB-K)/20)

where K is a calibration factor, which for NovaSAR-1 is -83 dB.

Intensity (or backscatter power) to/from Amplitude:

intensity = amplitude * amplitude
amplitude = sqrt(intensity)

Additional reference: https://forum.step.esa.int/t/what-stage-of-processing-requires-the-linear-to-from-db-command

Set up¶

Imports¶

In [1]:
# Common imports and settings
import os, sys, re
from pathlib import Path
from IPython.display import Markdown
import pandas as pd
pd.set_option("display.max_rows", None)
import xarray as xr
import numpy as np

# Datacube
import datacube
from datacube.utils.aws import configure_s3_access
import odc.geo.xr                             # https://github.com/opendatacube/odc-geo
from datacube.utils import masking            # https://github.com/opendatacube/datacube-core/blob/develop/datacube/utils/masking.py
from dea_tools.plotting import display_map    # https://github.com/GeoscienceAustralia/dea-notebooks/tree/develop/Tools

# Dask
import dask
from dask.distributed import Client, LocalCluster

# Basic plots
%matplotlib inline
# import matplotlib.pyplot as plt
# plt.rcParams['figure.figsize'] = [12, 8]

# Holoviews
# https://holoviz.org/tutorial/Composing_Plots.html
# https://holoviews.org/user_guide/Composing_Elements.html
import hvplot.xarray
import panel as pn
from datashader import reductions
from holoviews import opts
In [2]:
# EASI defaults
# These are convenience functions so that the notebooks in this repository work in all EASI deployments (comment-out if not required)

repo = Path.home() / 'easi-notebooks'    # Change i s
if str(repo) not in sys.path:
    sys.path.append(str(repo))

# from easi_tools import EasiDefaults  # Convenience
from easi_tools import initialize_dask, xarray_object_size, heading  # Useful

EASI defaults¶

These default values are configured for each EASI instance. They help us to use the same training notebooks in each EASI instance. You may find some of the functions convenient for your work or you can easily override the values in your copy of this notebook.

In [3]:
# TODO update for a NovaSAR product or set of NovaSAR products

# easi = EasiDefaults()

# family = 'sentinel-1'
# product = easi.product(family)   # 'sentinel1_grd_gamma0_20m'
# display(Markdown(f'Default {family} product for "{easi.name}": [{product}]({easi.explorer}/products/{product})'))

Dask cluster¶

Using a local Dask cluster is a good habit to get into. It can simplify loading and processing of data in many cases, and it provides a dashboard that shows the loading/processing progress.

To learn more about Dask see the set of dask notebooks.

In [4]:
# Local cluster
# cluster, client = initialize_dask(workers=4)  # EasiDefaults not configured for dev
# display(client)

# Damien's dask dashboard user path not available in dev (yet)
cluster = LocalCluster(n_workers=4)
client = Client(cluster)
server = f'https://hub.dev.easi-eo.solutions'  # Or replace if not using EasiDefaults
user = os.environ.get('JUPYTERHUB_SERVICE_PREFIX')  # Current user
dask.config.set({"distributed.dashboard.link": f'{server}{user}' + "proxy/{port}/status"})  # port is evaluated by dask
display(client)

# Or use Dask Gateway - this may take a few minutes
# cluster, client = initialize_dask(use_gateway=True, workers=4)
# display(client)

Client

Client-7227cb23-3ab9-11f0-8cd4-31d68f33f52d

Connection method: Cluster object Cluster type: distributed.LocalCluster
Dashboard: https://hub.dev.easi-eo.solutions/user/matt.paget/proxy/8787/status

Cluster Info

LocalCluster

c920389b

Dashboard: https://hub.dev.easi-eo.solutions/user/matt.paget/proxy/8787/status Workers: 4
Total threads: 4 Total memory: 8.00 GiB
Status: running Using processes: True

Scheduler Info

Scheduler

Scheduler-a5fc89de-ccaf-4b78-bc23-db59f372d98e

Comm: tcp://127.0.0.1:40507 Workers: 4
Dashboard: https://hub.dev.easi-eo.solutions/user/matt.paget/proxy/8787/status Total threads: 4
Started: Just now Total memory: 8.00 GiB

Workers

Worker: 0

Comm: tcp://127.0.0.1:38439 Total threads: 1
Dashboard: https://hub.dev.easi-eo.solutions/user/matt.paget/proxy/39089/status Memory: 2.00 GiB
Nanny: tcp://127.0.0.1:46149
Local directory: /tmp/dask-scratch-space/worker-0cc9n8eb

Worker: 1

Comm: tcp://127.0.0.1:39667 Total threads: 1
Dashboard: https://hub.dev.easi-eo.solutions/user/matt.paget/proxy/34485/status Memory: 2.00 GiB
Nanny: tcp://127.0.0.1:42805
Local directory: /tmp/dask-scratch-space/worker-qwu4jt6i

Worker: 2

Comm: tcp://127.0.0.1:45165 Total threads: 1
Dashboard: https://hub.dev.easi-eo.solutions/user/matt.paget/proxy/41065/status Memory: 2.00 GiB
Nanny: tcp://127.0.0.1:36703
Local directory: /tmp/dask-scratch-space/worker-2w9mf8oi

Worker: 3

Comm: tcp://127.0.0.1:39371 Total threads: 1
Dashboard: https://hub.dev.easi-eo.solutions/user/matt.paget/proxy/45181/status Memory: 2.00 GiB
Nanny: tcp://127.0.0.1:41163
Local directory: /tmp/dask-scratch-space/worker-oble24_c

ODC database¶

Connect to the ODC database. Configure the environment and low-level tools to read from AWS buckets.

In [5]:
dc = datacube.Datacube()

# Access AWS "requester-pays" buckets
# This is necessary for reading data from most third-party AWS S3 buckets such as for Landsat and Sentinel-2
configure_s3_access(aws_unsigned=False, requester_pays=True, client=client);
/env/lib/python3.12/site-packages/datacube/cfg/api.py:301: UserWarning: The 'default_environment' setting in the 'user' section is no longer supported - please refer to the documentation for more information
  warnings.warn(
/env/lib/python3.12/site-packages/datacube/cfg/opt.py:106: ODC2DeprecationWarning: Config being passed in by legacy environment variable $DB_USERNAME. Please use $ODC_DEFAULT_DB_USERNAME instead.
  warnings.warn(
/env/lib/python3.12/site-packages/datacube/cfg/opt.py:106: ODC2DeprecationWarning: Config being passed in by legacy environment variable $DB_PASSWORD. Please use $ODC_DEFAULT_DB_PASSWORD instead.
  warnings.warn(
/env/lib/python3.12/site-packages/datacube/cfg/opt.py:106: ODC2DeprecationWarning: Config being passed in by legacy environment variable $DB_HOSTNAME. Please use $ODC_DEFAULT_DB_HOSTNAME instead.
  warnings.warn(
/env/lib/python3.12/site-packages/datacube/cfg/opt.py:106: ODC2DeprecationWarning: Config being passed in by legacy environment variable $DB_PORT. Please use $ODC_DEFAULT_DB_PORT instead.
  warnings.warn(
/env/lib/python3.12/site-packages/datacube/cfg/opt.py:106: ODC2DeprecationWarning: Config being passed in by legacy environment variable $DB_DATABASE. Please use $ODC_DEFAULT_DB_DATABASE instead.
  warnings.warn(

Example query¶

Change any of the parameters in the query object below to adjust the location, time, projection, or spatial resolution of the returned datasets.

Use the Explorer interface to check the temporal and spatial coverage for each product.

In [6]:
# Explorer link
# display(Markdown(f'See: {easi.explorer}/products/{product}'))

# # EASI defaults
# display(Markdown(f'#### Location: {easi.location}'))
# latitude_range = easi.latitude
# longitude_range = easi.longitude
# time_range = easi.time

# Or set your own latitude / longitude
# Australia GWW
# latitude_range = (-33, -32.6)
# longitude_range = (120.5, 121)
# time_range = ('2020-01-01', '2020-01-31')

# Mt Gambier
# https://explorer.dev.easi-eo.solutions/products/novasar_ceosard_scd_vv_hh_hv/datasets/a522f17a-8537-52e1-b574-d9ef290c151f
# https://novasar.dev.easi-eo.solutions/collections/novasar_ceosard_scd_vv_hh_hv/items/a522f17a-8537-52e1-b574-d9ef290c151f
product = 'novasar_ceosard_scd_vv_hh_hv'
longitude_range = (140.2, 141.0)
latitude_range = (-38.1, -37.5)
time_range = '2021-06-12'

products = [
    'novasar_ceosard_grd_hh',  # NovaSAR-1 CEOS-ARD GRD HH polarization
    'novasar_ceosard_grd_vv',  # NovaSAR-1 CEOS-ARD GRD VV polarization
    'novasar_ceosard_scd_hh',  # NovaSAR-1 CEOS-ARD SCD HH polarization
    'novasar_ceosard_scd_hh_hv',  # NovaSAR-1 CEOS-ARD SCD HH,HV polarizations
    'novasar_ceosard_scd_hv',  # NovaSAR-1 CEOS-ARD SCD HV polarization
    'novasar_ceosard_scd_vv',  # NovaSAR-1 CEOS-ARD SCD VV polarization
    'novasar_ceosard_scd_vv_hh',  # NovaSAR-1 CEOS-ARD SCD VV,HH polarizations
    'novasar_ceosard_scd_vv_hh_hv',  # NovaSAR-1 CEOS-ARD SCD VV,HH,HV polarizations
]

query = {
    'product': 'novasar_ceosard_scd_vv_hh_hv',       # Product name
    'measurements': ['vv', 'hh', 'hv', 'angle', 'scatteringarea', 'gammatosigmaratio', 'mask'],  # All bands for testing but using alias labels for consistency with S1
    'x': longitude_range,     # "x" axis bounds
    'y': latitude_range,      # "y" axis bounds
    'time': time_range,       # Any parsable date strings
}

# Convenience function to display the selected area of interest
display_map(longitude_range, latitude_range)
Out[6]:
Make this Notebook Trusted to load map: File -> Trust Notebook

Load data¶

In [7]:
# Target xarray parameters
# - Select a set of measurements to load
# - output CRS and resolution
# - Usually we group input scenes on the same day to a single time layer (groupby)
# - Select a reasonable Dask chunk size (this should be adjusted depending on the
#   spatial and resolution parameters you choose
load_params = {
    'group_by': 'solar_day',                        # Scene grouping
    'dask_chunks': {'latitude':2048, 'longitude':2048},      # Dask chunks
}

# Load data
data = dc.load(**(query | load_params))
display(xarray_object_size(data))
display(data)
'Dataset size: 217.62 MB'
<xarray.Dataset> Size: 228MB
Dimensions:            (time: 1, latitude: 3001, longitude: 4001)
Coordinates:
  * time               (time) datetime64[ns] 8B 2021-06-12T01:26:39.429245
  * latitude           (latitude) float64 24kB -37.5 -37.5 -37.5 ... -38.1 -38.1
  * longitude          (longitude) float64 32kB 140.2 140.2 ... 141.0 141.0
    spatial_ref        int32 4B 4326
Data variables:
    vv                 (time, latitude, longitude) uint16 24MB dask.array<chunksize=(1, 2048, 2048), meta=np.ndarray>
    hh                 (time, latitude, longitude) uint16 24MB dask.array<chunksize=(1, 2048, 2048), meta=np.ndarray>
    hv                 (time, latitude, longitude) uint16 24MB dask.array<chunksize=(1, 2048, 2048), meta=np.ndarray>
    angle              (time, latitude, longitude) float32 48MB dask.array<chunksize=(1, 2048, 2048), meta=np.ndarray>
    scatteringarea     (time, latitude, longitude) float32 48MB dask.array<chunksize=(1, 2048, 2048), meta=np.ndarray>
    gammatosigmaratio  (time, latitude, longitude) float32 48MB dask.array<chunksize=(1, 2048, 2048), meta=np.ndarray>
    mask               (time, latitude, longitude) uint8 12MB dask.array<chunksize=(1, 2048, 2048), meta=np.ndarray>
Attributes:
    crs:           EPSG:4326
    grid_mapping:  spatial_ref
xarray.Dataset
    • time: 1
    • latitude: 3001
    • longitude: 4001
    • time
      (time)
      datetime64[ns]
      2021-06-12T01:26:39.429245
      units :
      seconds since 1970-01-01 00:00:00
      array(['2021-06-12T01:26:39.429245000'], dtype='datetime64[ns]')
    • latitude
      (latitude)
      float64
      -37.5 -37.5 -37.5 ... -38.1 -38.1
      units :
      degrees_north
      resolution :
      -0.0002
      crs :
      EPSG:4326
      array([-37.5   , -37.5002, -37.5004, ..., -38.0996, -38.0998, -38.1   ])
    • longitude
      (longitude)
      float64
      140.2 140.2 140.2 ... 141.0 141.0
      units :
      degrees_east
      resolution :
      0.0002
      crs :
      EPSG:4326
      array([140.2   , 140.2002, 140.2004, ..., 140.9996, 140.9998, 141.    ])
    • spatial_ref
      ()
      int32
      4326
      spatial_ref :
      GEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 (Transit)"],MEMBER["World Geodetic System 1984 (G730)"],MEMBER["World Geodetic System 1984 (G873)"],MEMBER["World Geodetic System 1984 (G1150)"],MEMBER["World Geodetic System 1984 (G1674)"],MEMBER["World Geodetic System 1984 (G1762)"],MEMBER["World Geodetic System 1984 (G2139)"],MEMBER["World Geodetic System 1984 (G2296)"],ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],USAGE[SCOPE["Horizontal component of 3D system."],AREA["World."],BBOX[-90,-180,90,180]],ID["EPSG",4326]]
      crs_wkt :
      GEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 (Transit)"],MEMBER["World Geodetic System 1984 (G730)"],MEMBER["World Geodetic System 1984 (G873)"],MEMBER["World Geodetic System 1984 (G1150)"],MEMBER["World Geodetic System 1984 (G1674)"],MEMBER["World Geodetic System 1984 (G1762)"],MEMBER["World Geodetic System 1984 (G2139)"],MEMBER["World Geodetic System 1984 (G2296)"],ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],USAGE[SCOPE["Horizontal component of 3D system."],AREA["World."],BBOX[-90,-180,90,180]],ID["EPSG",4326]]
      semi_major_axis :
      6378137.0
      semi_minor_axis :
      6356752.314245179
      inverse_flattening :
      298.257223563
      reference_ellipsoid_name :
      WGS 84
      longitude_of_prime_meridian :
      0.0
      prime_meridian_name :
      Greenwich
      geographic_crs_name :
      WGS 84
      horizontal_datum_name :
      World Geodetic System 1984 ensemble
      grid_mapping_name :
      latitude_longitude
      GeoTransform :
      140.199900000000013733369997 0.000200000000000000009584 0 -37.499899999999996680344339 0 -0.000200000000000000009584
      array(4326, dtype=int32)
    • vv
      (time, latitude, longitude)
      uint16
      dask.array<chunksize=(1, 2048, 2048), meta=np.ndarray>
      units :
      amplitude
      nodata :
      0
      add_offset :
      0
      scale_factor :
      1
      crs :
      EPSG:4326
      grid_mapping :
      spatial_ref
      Array Chunk
      Bytes 22.90 MiB 8.00 MiB
      Shape (1, 3001, 4001) (1, 2048, 2048)
      Dask graph 4 chunks in 1 graph layer
      Data type uint16 numpy.ndarray
      4001 3001 1
    • hh
      (time, latitude, longitude)
      uint16
      dask.array<chunksize=(1, 2048, 2048), meta=np.ndarray>
      units :
      amplitude
      nodata :
      0
      add_offset :
      0
      scale_factor :
      1
      crs :
      EPSG:4326
      grid_mapping :
      spatial_ref
      Array Chunk
      Bytes 22.90 MiB 8.00 MiB
      Shape (1, 3001, 4001) (1, 2048, 2048)
      Dask graph 4 chunks in 1 graph layer
      Data type uint16 numpy.ndarray
      4001 3001 1
    • hv
      (time, latitude, longitude)
      uint16
      dask.array<chunksize=(1, 2048, 2048), meta=np.ndarray>
      units :
      amplitude
      nodata :
      0
      add_offset :
      0
      scale_factor :
      1
      crs :
      EPSG:4326
      grid_mapping :
      spatial_ref
      Array Chunk
      Bytes 22.90 MiB 8.00 MiB
      Shape (1, 3001, 4001) (1, 2048, 2048)
      Dask graph 4 chunks in 1 graph layer
      Data type uint16 numpy.ndarray
      4001 3001 1
    • angle
      (time, latitude, longitude)
      float32
      dask.array<chunksize=(1, 2048, 2048), meta=np.ndarray>
      units :
      deg
      nodata :
      0
      crs :
      EPSG:4326
      grid_mapping :
      spatial_ref
      Array Chunk
      Bytes 45.80 MiB 16.00 MiB
      Shape (1, 3001, 4001) (1, 2048, 2048)
      Dask graph 4 chunks in 1 graph layer
      Data type float32 numpy.ndarray
      4001 3001 1
    • scatteringarea
      (time, latitude, longitude)
      float32
      dask.array<chunksize=(1, 2048, 2048), meta=np.ndarray>
      units :
      m^2
      nodata :
      0
      crs :
      EPSG:4326
      grid_mapping :
      spatial_ref
      Array Chunk
      Bytes 45.80 MiB 16.00 MiB
      Shape (1, 3001, 4001) (1, 2048, 2048)
      Dask graph 4 chunks in 1 graph layer
      Data type float32 numpy.ndarray
      4001 3001 1
    • gammatosigmaratio
      (time, latitude, longitude)
      float32
      dask.array<chunksize=(1, 2048, 2048), meta=np.ndarray>
      units :
      1
      nodata :
      0
      crs :
      EPSG:4326
      grid_mapping :
      spatial_ref
      Array Chunk
      Bytes 45.80 MiB 16.00 MiB
      Shape (1, 3001, 4001) (1, 2048, 2048)
      Dask graph 4 chunks in 1 graph layer
      Data type float32 numpy.ndarray
      4001 3001 1
    • mask
      (time, latitude, longitude)
      uint8
      dask.array<chunksize=(1, 2048, 2048), meta=np.ndarray>
      units :
      bits
      nodata :
      0
      flags_definition :
      {'mask': {'bits': [0, 1, 2, 3, 4, 5, 6, 7], 'values': {'0': 'InvalidData', '1': 'ValidData', '5': 'Layover', '17': 'Shadow'}, 'description': 'Data Mask Image'}}
      crs :
      EPSG:4326
      grid_mapping :
      spatial_ref
      Array Chunk
      Bytes 11.45 MiB 4.00 MiB
      Shape (1, 3001, 4001) (1, 2048, 2048)
      Dask graph 4 chunks in 1 graph layer
      Data type uint8 numpy.ndarray
      4001 3001 1
    • time
      PandasIndex
      PandasIndex(DatetimeIndex(['2021-06-12 01:26:39.429245'], dtype='datetime64[ns]', name='time', freq=None))
    • latitude
      PandasIndex
      PandasIndex(Index([   -37.5, -37.5002, -37.5004, -37.5006, -37.5008,  -37.501, -37.5012,
             -37.5014, -37.5016, -37.5018,
             ...
             -38.0982, -38.0984, -38.0986, -38.0988,  -38.099, -38.0992, -38.0994,
             -38.0996, -38.0998,    -38.1],
            dtype='float64', name='latitude', length=3001))
    • longitude
      PandasIndex
      PandasIndex(Index([140.20000000000002, 140.20020000000002, 140.20040000000003,
                       140.2006, 140.20080000000002, 140.20100000000002,
             140.20120000000003,           140.2014,           140.2016,
             140.20180000000002,
             ...
             140.99820000000003,           140.9984,           140.9986,
             140.99880000000002, 140.99900000000002, 140.99920000000003,
                       140.9994, 140.99960000000002, 140.99980000000002,
             141.00000000000003],
            dtype='float64', name='longitude', length=4001))
  • crs :
    EPSG:4326
    grid_mapping :
    spatial_ref
In [8]:
# When happy with the shape and size of chunks, persist() the result
data = data.persist()

Conversion and helper functions¶

In [9]:
# These functions use numpy, which should be satisfactory for most notebooks.
# Calculations for larger or more complex arrays may require Xarray's "ufunc" capability.
# https://docs.xarray.dev/en/stable/examples/apply_ufunc_vectorize_1d.html
#
# Apply numpy.log10 to the DataArray
# log10_data = xr.apply_ufunc(np.log10, data)

def dn_to_decibel(da: 'xr.DataArray', K=0):
    """Return an array converted to dB values"""
    xx = da.where(da > 0, np.nan)  # Set values <= 0 to NaN
    xx = 20*np.log10(xx) + K
    xx.attrs.update({"units": "dB"})
    return xx

def decibel_to_dn(da: 'xr.DataArray', K=0):
    """Return an array converted to intensity values"""
    xx = np.power(10, (da-K)/20.0)
    xx.attrs.update({"units": "intensity"})
    return xx

def dn_to_amplitude(da: 'xr.DataArray'):
    """Return an array converted to amplitude values"""
    # TODO get scaling-factor from measurement (ODC database)
    xx = da / 14125.3754
    xx.attrs.update({"units": "amplitude"})
    return xx

def select_valid_time_layers(ds: 'xarray', percent: float = 5):
    """Select time layers that have at least a given percentage of valid data (e.g., >=5%)

    Example usage:
      selected = select_valid_time_layers(ds, percent=5)
      filtered == ds.sel(time=selected)
    """
    spatial_dims = ds.odc.spatial_dims
    return ds.count(dim=spatial_dims).values / (ds.sizes[spatial_dims[0]]*ds.sizes[spatial_dims[1]]) >= (percent/100.0)

# Examples to check that the intensity to/from dB functions work as expected
# xx = data.vv.isel(time=0,latitude=np.arange(0, 5),longitude=np.arange(0, 5))
# xx[0] = 0
# xx[1] = -0.001
# display(xx.values)
# yy = dn_to_decibel(xx)
# display(yy.values)
# zz = decibel_to_dn(yy)
# display(zz.values)
In [10]:
# hvPlot convenience functions
def make_image(ds: 'xarray', frame_height=300, **kwargs):
    """Return a Holoviews DynamicMap (image) object that can be displayed or combined"""
    spatial_dims = ds.odc.spatial_dims
    defaults = dict(
        cmap="Greys_r",
        y = spatial_dims[0], x = spatial_dims[1],
        groupby = 'time',
        rasterize = True,
        geo = True,
        robust = True,
        frame_height = frame_height,
        clabel = ds.attrs.get('units', None),
    )
    defaults.update(**kwargs)
    return ds.hvplot.image(**defaults)

def rgb_image(ds: 'xarray', frame_height=300, **kwargs):
    """Return a Holoviews DynamicMap (RBG image) object that can be displayed or combined"""
    spatial_dims = ds.odc.spatial_dims
    defaults = dict(
        bands='band',
        y = spatial_dims[0], x = spatial_dims[1],
        groupby = 'time',
        rasterize = True,
        geo = True,
        robust = True,
        frame_height = frame_height,
    )
    defaults.update(**kwargs)
    return ds.hvplot.rgb(**defaults)


from bokeh.models.tickers import FixedTicker
def mask_image(ds: 'xarray', frame_height=300, **kwargs):
    """Return a Holoviews DynamicMap (RBG image) object that can be displayed or combined"""
    # NovaSAR ARD mask
    color_def = [
        (0,  '#ff0004', 'InvalidData'),   # red
        (1,  '#eeeeee', 'ValidData'),   # light grey
        (5, '#ff52ff', 'Layover'),   # cyan
        (17,  '#774c0b', 'Shadow'),   # brown
    ]
    color_val = [x[0] for x in color_def]
    color_hex = [x[1] for x in color_def]
    color_txt = [f'{x[0]:2d}: {x[2]}' for x in color_def]
    color_lim = (min(color_val), max(color_val) + 1)
    bin_edges = color_val + [max(color_val) + 1]
    # bin_range = (color_val[0] + 0.5, color_val[-1] + 0.5)  # No idea why (0.5,11.5) works and (0,11) or (0,12) do not

    # Image options
    defaults = {
        'aggregator': reductions.mode(),
        'cmap': color_hex,
        'clim': color_lim,
        'colorbar': True,
        'frame_height': frame_height,
    }
    # Colorbar options for categories
    extra_opts = {
        'color_levels': bin_edges,
        'colorbar_opts': {
            'major_label_overrides': dict(zip(color_val, color_txt)),
            'major_label_text_align': 'left',
            'ticker': FixedTicker(ticks=color_val),
        },
    }
    defaults.update(**kwargs)
    return make_image(ds, **defaults).options(opts.Image(**extra_opts))  #.hist(bins=bin_edges)
In [11]:
# Optional time layer filter

# selected = select_valid_time_layers(data.vv, 10)  # Exclude time layers with less than 10% valid data
# data = data.sel(time=selected).persist()
In [12]:
# Add dB and amplitude values to the dataset

data['vv_db'] = dn_to_decibel(data.vv, K=-83).persist()
data['hh_db'] = dn_to_decibel(data.hh, K=-83).persist()
data['hv_db'] = dn_to_decibel(data.hv, K=-83).persist()

data['vv_amp'] = dn_to_amplitude(data.vv).persist()
data['hh_amp'] = dn_to_amplitude(data.hh).persist()
data['hv_amp'] = dn_to_amplitude(data.hv).persist()

data.mask.persist()
Out[12]:
<xarray.DataArray 'mask' (time: 1, latitude: 3001, longitude: 4001)> Size: 12MB
dask.array<dc_load_mask, shape=(1, 3001, 4001), dtype=uint8, chunksize=(1, 2048, 2048), chunktype=numpy.ndarray>
Coordinates:
  * time         (time) datetime64[ns] 8B 2021-06-12T01:26:39.429245
  * latitude     (latitude) float64 24kB -37.5 -37.5 -37.5 ... -38.1 -38.1 -38.1
  * longitude    (longitude) float64 32kB 140.2 140.2 140.2 ... 141.0 141.0
    spatial_ref  int32 4B 4326
Attributes:
    units:             bits
    nodata:            0
    flags_definition:  {'mask': {'bits': [0, 1, 2, 3, 4, 5, 6, 7], 'values': ...
    crs:               EPSG:4326
    grid_mapping:      spatial_ref
xarray.DataArray
'mask'
  • time: 1
  • latitude: 3001
  • longitude: 4001
  • dask.array<chunksize=(1, 2048, 2048), meta=np.ndarray>
    Array Chunk
    Bytes 11.45 MiB 4.00 MiB
    Shape (1, 3001, 4001) (1, 2048, 2048)
    Dask graph 4 chunks in 1 graph layer
    Data type uint8 numpy.ndarray
    4001 3001 1
    • time
      (time)
      datetime64[ns]
      2021-06-12T01:26:39.429245
      units :
      seconds since 1970-01-01 00:00:00
      array(['2021-06-12T01:26:39.429245000'], dtype='datetime64[ns]')
    • latitude
      (latitude)
      float64
      -37.5 -37.5 -37.5 ... -38.1 -38.1
      units :
      degrees_north
      resolution :
      -0.0002
      crs :
      EPSG:4326
      array([-37.5   , -37.5002, -37.5004, ..., -38.0996, -38.0998, -38.1   ])
    • longitude
      (longitude)
      float64
      140.2 140.2 140.2 ... 141.0 141.0
      units :
      degrees_east
      resolution :
      0.0002
      crs :
      EPSG:4326
      array([140.2   , 140.2002, 140.2004, ..., 140.9996, 140.9998, 141.    ])
    • spatial_ref
      ()
      int32
      4326
      spatial_ref :
      GEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 (Transit)"],MEMBER["World Geodetic System 1984 (G730)"],MEMBER["World Geodetic System 1984 (G873)"],MEMBER["World Geodetic System 1984 (G1150)"],MEMBER["World Geodetic System 1984 (G1674)"],MEMBER["World Geodetic System 1984 (G1762)"],MEMBER["World Geodetic System 1984 (G2139)"],MEMBER["World Geodetic System 1984 (G2296)"],ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],USAGE[SCOPE["Horizontal component of 3D system."],AREA["World."],BBOX[-90,-180,90,180]],ID["EPSG",4326]]
      crs_wkt :
      GEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 (Transit)"],MEMBER["World Geodetic System 1984 (G730)"],MEMBER["World Geodetic System 1984 (G873)"],MEMBER["World Geodetic System 1984 (G1150)"],MEMBER["World Geodetic System 1984 (G1674)"],MEMBER["World Geodetic System 1984 (G1762)"],MEMBER["World Geodetic System 1984 (G2139)"],MEMBER["World Geodetic System 1984 (G2296)"],ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],USAGE[SCOPE["Horizontal component of 3D system."],AREA["World."],BBOX[-90,-180,90,180]],ID["EPSG",4326]]
      semi_major_axis :
      6378137.0
      semi_minor_axis :
      6356752.314245179
      inverse_flattening :
      298.257223563
      reference_ellipsoid_name :
      WGS 84
      longitude_of_prime_meridian :
      0.0
      prime_meridian_name :
      Greenwich
      geographic_crs_name :
      WGS 84
      horizontal_datum_name :
      World Geodetic System 1984 ensemble
      grid_mapping_name :
      latitude_longitude
      GeoTransform :
      140.199900000000013733369997 0.000200000000000000009584 0 -37.499899999999996680344339 0 -0.000200000000000000009584
      array(4326, dtype=int32)
    • time
      PandasIndex
      PandasIndex(DatetimeIndex(['2021-06-12 01:26:39.429245'], dtype='datetime64[ns]', name='time', freq=None))
    • latitude
      PandasIndex
      PandasIndex(Index([   -37.5, -37.5002, -37.5004, -37.5006, -37.5008,  -37.501, -37.5012,
             -37.5014, -37.5016, -37.5018,
             ...
             -38.0982, -38.0984, -38.0986, -38.0988,  -38.099, -38.0992, -38.0994,
             -38.0996, -38.0998,    -38.1],
            dtype='float64', name='latitude', length=3001))
    • longitude
      PandasIndex
      PandasIndex(Index([140.20000000000002, 140.20020000000002, 140.20040000000003,
                       140.2006, 140.20080000000002, 140.20100000000002,
             140.20120000000003,           140.2014,           140.2016,
             140.20180000000002,
             ...
             140.99820000000003,           140.9984,           140.9986,
             140.99880000000002, 140.99900000000002, 140.99920000000003,
                       140.9994, 140.99960000000002, 140.99980000000002,
             141.00000000000003],
            dtype='float64', name='longitude', length=4001))
  • units :
    bits
    nodata :
    0
    flags_definition :
    {'mask': {'bits': [0, 1, 2, 3, 4, 5, 6, 7], 'values': {'0': 'InvalidData', '1': 'ValidData', '5': 'Layover', '17': 'Shadow'}, 'description': 'Data Mask Image'}}
    crs :
    EPSG:4326
    grid_mapping :
    spatial_ref

Plot the data¶

Note the different data ranges for plotting (clim) between vv, vh, intensity and dB.

  • Stronger co-polarisation (VV) indicates direct backscatter while stronger cross-polarisation (VH) may indicate a complex surface or volume scattering.
  • Intensity data are linear-scaled so can tend to disciminate across a range of backscatter returns.
  • Decibel data are log-scaled so can tend to discriminate high and low backscatter returns.
In [13]:
# VV, HH and HV (amplitude and dB) and Angle hvPlots

vv_plot = make_image(data.vv_amp, title='VV (amplitude)', clim=(0, 0.5))
hh_plot = make_image(data.hh_amp, title='HH (amplitude)', clim=(0, 0.5))
hv_plot = make_image(data.hv_amp, title='HV (amplitude)', clim=(0, 0.25))

vv_db_plot = make_image(data.vv_db, title='VV (dB)', clim=(-20, -5))
hh_db_plot = make_image(data.hh_db, title='HH (dB)', clim=(-20, -5))
hv_db_plot = make_image(data.hv_db, title='HV (dB)', clim=(-35, -10))

ia_plot = make_image(data.angle, title='Incidence angle')
sa_plot = make_image(data.scatteringarea, title='Scattering area')
gs_plot = make_image(data.gammatosigmaratio, title='Gamma to sigma ratio')
mask_plot = mask_image(data.mask, title='Mask')
In [14]:
# Arrange plots with linked axes and time slider. Adjust browser window width if required.

layout = pn.panel(
    (vv_plot + hh_plot + hv_plot + \
     vv_db_plot + hh_db_plot + hv_db_plot + \
     ia_plot + sa_plot + gs_plot + mask_plot).cols(3),
    widget_location='top',
)
print(layout)  # Helpful to see how the hvplot is constructed
layout
Column
    [0] WidgetBox(align=('center', 'start'))
        [0] Select(name='time (seconds s..., options={'2021-06-12 0...}, value=numpy.datetime64('2021-06-...)
    [1] HoloViews(Layout, widget_location='top')
Out[14]:

-------------¶

The remainder of this notebook is yet to be updated for NovaSAR

Plot a histogram of the dB data¶

A histogram can help separate water from land features. Here we show a histogram for the HV (db) channel for all time layers.

  • If the histogram shows two clear peaks then a value between the peaks could be used as a water / land threshold
  • If not then try selected time layers, a different area of interest, or other channels or combinations.
In [15]:
# vals, bins, hist_plot = data.hv_db.plot.hist(bins=np.arange(-30, 0, 1), color='red')  # Matplotlib
hist_plot = data.hv_db.hvplot.hist(bins=np.arange(-30, 0, 1), color='red', title='Combined times', height=400)  # hvPlot

print(hist_plot)  # Helpful to see how the hvplot is constructed
hist_plot
:NdOverlay   [Variable]
   :Histogram   [hv_db]   (Count)
Out[15]:

Make an RGB image¶

A common strategy to create an RGB colour composite image for SAR data from two channels is to use the ratio of the channels to represent the third colour. Here we choose

To create an RGB colour composite image we can use the ratio of VH and VV to represent a third channel. Here we choose

  • Red = VHV ... complex scattering
  • Green = VV ... direct scattering
  • Blue = VH/VV ... relatively more complex than direct
In [ ]:
# Add the vh/vv band to represent 'blue'
data['vh_vv'] = data.vh / data.vv

# Scale the measurements by their median so they have a similar range for visualization
spatial_dims = data.odc.spatial_dims
data['vh_scaled'] = data.vh / data.vh.median(dim=spatial_dims).persist()
data['vv_scaled'] = data.vv / data.vv.median(dim=spatial_dims).persist()
data['vh_vv_scaled'] = data.vh_vv / data.vh_vv.median(dim=spatial_dims).persist()

# odc-geo function
rgb_data = data.odc.to_rgba(bands=['vh_scaled','vv_scaled','vh_vv_scaled'], vmin=0, vmax=2)
In [ ]:
# As subplots
# rgb_plot = rgb_image(
#     rgb_data,
# ).layout().cols(4)

# As movie. Select "loop" and use "-" button to adjust the speed to allow for rendering. After a few cycles the images should play reasonably well.
rgb_plot = rgb_image(
    rgb_data,
    precompute = True,
    widget_type='scrubber', widget_location='bottom',
    frame_height = 500,
)

print(rgb_plot)  # Helpful to see how the hvplot is constructed
rgb_plot

Export to Geotiffs¶

Recall that to write a dask dataset to a file requires the dataset to be .compute()ed. This may result in a large memory increase on your JupyterLab node if the area of interest is large enough, which in turn may kill the kernel. If so then skip this step, choose a smaller area or find a different way to export data.

In [ ]:
# Make a directory to save outputs to
target = Path.home() / 'output'
if not target.exists(): target.mkdir()

def write_band(ds, varname):
    """Write the variable name of the xarray dataset to a Geotiff file for each time layer"""
    for i in range(len(ds.time)):
        date = ds[varname].isel(time=i).time.dt.strftime('%Y%m%d').data
        fname = f'{target}/example_sentinel-1_{varname}_{date}.tif'
        single = ds[varname].isel(time=i).compute()
        single.odc.write_cog(
            fname=fname,
            overwrite=True,
        )
        print(f'Wrote: {fname}')
        
write_band(data, 'vv')
write_band(data, 'vh')
In [ ]: